Spectral density of the Hubbard-model by the continued fraction method 
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We present the continued fraction method (CFM) as a new microscopic approximation to the 
spectral density of the Hubbard model in the correlated metal phase away from half filling. The 
quantity expanded as a continued fraction is the single particle Green function. Leading spectral 
moments are taken into account through a set of real expansion coefficients, as known from the 
projection technique. The new aspect is to add further stages to the continued fraction, with complex 
coefficients, thus defining a terminator function. This enables us to treat the entire spectral range of 
the Green function on equal footing and determine the energy scale of the Fermi liquid quasiparticles 
by minimizing the total energy. The solution is free of phenomenological parameters and remains 
well defined in the strong coupling limit, near the doping controlled metal-insulator transition. Our 
results for the density of states agree reasonably with several variants of the dynamical mean field 
theory. The CFM requires minimal numerical effort and can be generalized in several ways that are 
interesting for applications to real materials. 

PACS numbers: 71.10.Fd, 75.20.Hr 
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I. INTRODUCTION 

The Hubbard Hamiltoniani is certainly the most im- 
portant model in the field of strongly correlated electrons. 
The spectral function for the addition or removal of a sin- 
gle electron near half filling serves as a paradigm for the 
excitation spectrum of highly correlated electrons in the 
vicinity of a Mott transition. It was a great success of 
the dynamical mean field theory (DMFT) to connect 
the high- and low-energy parts of the spectral function 
in a non-perturbative solution for arbitrary interaction 
strength. Especially, it was confirmed that the coherent 
low energy excitations in the metallic phase follow the 
same dynamics as the Kondo resonance in the Anderson 
impurity-model, the other generic Hamiltonian for corre- 
lated electrons that is much better understood^ 

The analog to the Kondo resonance in the impurity 
model is a quasiparticle (QP) band in the lattice model. 
Both straddle the chemical potential fx, i.e. the lowest 
excitations are gapless. In the strong coupling limit, the 
spectral weight Z of the QP band is small, relative to 
two sidebands, further removed from fx. These sidebands 
are called the Hubbard bands because they are roughly 
reminiscent of the Hubbard-I solution. 1 In the doping 
controlled regime^ one sideband always overlaps with the 
QPs, the other represents true high energy excitations 
across the correlation gap. 

Hubbard-I is an approximation close to the atomic 
limit, but nevertheless taking exact spectral moments of 
the itinerant propagator up to the second order into ac- 
count. It can be considered the ancestor of the projec- 
tion technique^ which systematically incorporates spec- 
tral moments of higher order. These approximations have 
severe deficiencies in the low energy sector, unless the mo- 



ment series can be effectively summed up. In particular, 
generating a third pole in the spectral function from the 
high energy side alone leads to uncontrolled results. 

When the density of states (DOS), as obtained within 
DMFT, is resolved with respect to the wavenumber k, 
more details about the coexistence of this Kondo reso- 
nance with atomic like features in a lattice system are 
revealed. The lowest excitations are true Fermi liquid 
(FL) QP's: (i) The finite DOS at e = fj, corresponds to 
longlived excitations, (ii) These are located in fc-space on 
a Fermi surface (FS) that satisfies Luttinger's theorem^ 
(iii) As function of the distance k — Uf from the FS, the 
excitation energy has a linear, strongly reduced disper- 
sion, (iv) The damping is quadratic in k—kp but strongly 
enhanced, meaning that the linear and quadratic term 
are of the same order at a very small energy scale, the 
coherence energy A*. The two atomic like excitations 
turn out to be strongly damped, even when k is on the 
FS. Their peaks disperse with k but spectral tails spread 
over the entire bandwidth. 

The DMFT thus unites atomic and itinerant features 
in a non perturbative approximation. It is exact only 
in dimension d = oo. As a generic scenario, it is ex- 
pected to hold down to d = 2, albeit with the caveat 
that the DMFT suppresses additional structure due to 
bosonic couplings. Earlier approximations at finite d al- 
ready yielded QP's^& and established a connection to 
the Kondo effect. 3 The high prestige of the DMFT is 
due to its ability to produce a selfconsistent, numerically 
manageable approximation to the spectral function for all 
energies, in particular to the parameter Z that governs 
the low energy sector. This has opened a path to realis- 
tic modeling of correlated materials beyond the Hubbard 
model in such methods as LDA+DMFT. 9 
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It is nevertheless desirable for several reasons to pursue 
alternative methods in parallel. Firstly, a fc-independent 
selfenergy, which is the proper result at d = oo, does not 
allow to explain phenomena that depend on the different 
symmetry directions in the Brillouin zone, especially in 
high temperature superconductors and other low dimen- 
sional systems. Cluster extensions of the DMFT go in 
the direction of lifting this restriction^ but these gen- 
eralizations are numerically even more demanding than 
the DMFT itself. The precise solution of a manybody 
Kondo problem is required at each iteration step towards 
selfconsistency. In practice, when designing the "impu- 
rity solver" , a trade-off exists between improving the low 
energy, low temperature solution and exactly satisfying 
global sumrules. Such numerical problems are presently 
a bottleneck for extensions of the DMFT to larger clus- 
ters or to LDA+DMFT with charge transfer into lig- 
and bands. A variational aspect was recently found, 
which may allow to circumvent some of the numerical 
problems^- 

In this paper, we present a continued fraction method 
(CFM) and implement it for the doping controlled metal- 
lic regime near the Mott transition. Similar to other re- 
cent attempts (i24i2iiiiiSii& we start from the projection 
technique, applied to the fc-resolved single particle Green 
function (GF). The notations are introduced in section 
ITT1 In section ITTT1 the connections between the moment- 
and continued fraction- (CF) expansion as well as the 
Pade approximant (PA) are established. In the PA, qual- 
itatively important features of the macroscopic system, 
such as damping, are missing. They can only be cap- 
tured by resummation of the CF to infinite order. The 
concept of a terminator function (TF), by which an ap- 
proximate resummation is achieved, is common to many 
methods based on the CF. As a general scheme, we de- 
fine our CFM by allowing only such TF's that preserve 
the structure of a truncated CF, however with complex 
coefficients. Useful recursion relations, that are proper- 
ties of PAs, can thus be carried over and the solution for 
the GF can be constrained by high as well as low energy 
sumrules. 

Previous solutions obtained with this ansatai^ were 
partly phenomenological, because the strong coupling 
renormalization Z needed to be inferred from a separate 
Gutzwiller approximation, or else was left open for fitting 
to experiments ! 17 i 18 A closed solution is now achieved by 
minimizing the total energy in the presence of sumrules, 
for which the necessary selfconsistency loops are intro- 
duced. The selfconsistent Z falls below the Gutzwiller 
value and has a doping dependence close to that for the 
exact Kondo scaled This is now a true microscopic ap- 
proximation, depending only on the parameters in the 
Hamiltonian. 

In sections llVl Ivl and I VII we have investigated the TF's 
that correspond to adding one or two stages with com- 
plex coefficients to the CF. We show how the FS singular- 
ity, the enclosed Luttinger volume in /c-space and the FL 
damping can be modeled rigorously. We assess the qual- 



ity of our approximations by comparing the DOS with 
the DMFT result for two variants of the impurity solver, 
namely the numerical renormalization group (NRG) 20 
and the non-crossing approximation (NCA)iSi 

The success of the CFM with respect to the Hubbard 
model allows to draw some optimistic conclusions about 
possible generalizations towards more realistic models, 
describing correlation effects in a multiband electronic 
environment. This will be outlined as part of the conclu- 
sions. 



II. HAMILTONIAN, GREEN FUNCTION AND 
GENERALITIES ABOUT CONTINUED 
FRACTIONS 

The Hubbard model for a grand canonical ensemble of 
electrons on a lattice of N sites (N — > oo) is written in 
the usual notations 

H = H-(iN = ^(£ fe - ii)c\ a c ka + u J2 mrmi ■ CO 

k.a i 

The kinetic energy consists of itinerant Bloch states with 
energies E k and wavenumber k, running through one 
Brillouin zone. The bandwidth is 2D and, when not 
specified otherwise, D is used as unit of both energy and 
frequency (h= 1). We formulate the method for an arbi- 
trary density of Bloch states. Numerical examples later 
on will be calculated for a semi-elliptic density. 

The chemical potential fx is selfconsistently determined 
to satisfy the condition 



n={N)/N = 2m = ^2(4a^), (2) 

cr 

where the filling factor n (0 < n < 2) is part of the input. 
The chemical potential for the U = limit is designated 
as /io- For U ^ 0, the right hand side is calculated with 
our method. The overline and the bracket signify Bril- 
louin zone average and ensemble average, respectively. 
The filling factor per spin direction in the spin degener- 
ate phase is m = n/2. 

We approximate the advanced single particle GF 

G(k,w)=i [ e^([c k<7 (t),cl(0)}}, (3) 

from which the momentum distribution (ct Cfc CT ) and 
other observables are calculated. The spin index is 
dropped in the unpolarized phase. The time dependent 
fermionic destruction operator Cka(t) is in the Hciscn- 
berg representation with H and the square bracket is the 
anticommutator. The complex frequency u> has ji as ori- 
gin. Asymptotically, for large u>, we have G{k,uj) ~ l/u>. 
The coefficient 1 reflects the moment M = 1 or spectral 
norm, as required by the Pauli principle. For this relation 
between the leading coefficient and the norm to remain 
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valid in an approximation, it is necessary and sufficient 
to conserve the Herglotz property. In the case of the 
advanced GF, it means that the relation lmG(k,u>) > 
must be obeyed throughout the entire halfplanc Imuj < 0. 
The physical meaning of the Herglotz property is causal- 
ity and it automatically entails the existence of Kramers- 
Kronig relations between the real and imaginary parts. A 
great advantage of our method is the possibility to make 
straightforward evaluations along the real axis. Since 
this limit has to be approached from within the domain 
of analyticity, the notation uj — e — i0 + with real energy 
e = E — /.t is introduced. The A-resolvcd spectral function 
is 



A(k,e) = -ImG(fc,e-iO+). 



(4) 



At U = it has a single sharp peak at the excitation 
energy 



Vk 



E k - (j, cc (k - k F ), 



(5) 



which also serves to measure distance in fc-space, at least 
in the vivinity of the FS. 

The CF expansion, on which our method is based 
in a crucial way, has already a long tradition in solid 
state physics, in the one electron problem with disorder 23 
as well as in the many electron problem^ The CF is 
generated by various procedures like tridiagonalization, 
recursion- or Lanczos-methods. The Hubbard-I GF is the 
simplest example of a CF that has been truncated at low 
order. The exact GF for the Hubbard model on finite 
clusters is a CF which naturally ends at very high order. 
The CF for the infinite system does not end. Proper- 
ties of the thermodynamic limit, such as damping due to 
electron-electron scattering, emerge only after resumma- 
tion of the CF. Approximate resummation is achieved by 
the TF, an analytic function which also has the Herglotz 
property. 

A well chosen TF is thus expected to bring two im- 
provements to the approximation for the GF in dimen- 
sions d > 2: (i) From a set of discrete, more or less intense 
and more or less densely spaced Dirac peaks emerges the 
final shape of the continuous spectral density (see Ref . |2{| 
for tight-binding like models and Ref. 26 for strongly cor- 
related electrons), (ii) A Fermi surface (FS) discontinu- 
ity emerges in the momentum distribution {c\ a Ck a ) at 
temperatures below the strong coupling energy scale A*. 
The FL discontinuity and the correct FS volume will be 
incorporated in our ansatz. This means, we take the Lut- 
tinger theorem for granted and use it as a principle, even 
for strong coupling where there is no rigorous proof. The 
energy A* then comes out as part of the selfconsistent 
solution. 



III. HIGH-ENERGY PART 

The first moment or center of gravity of G(k, uj) is 

u>i = Ek + mil — /i. (6) 



It disperses like the unrenormalized Bloch energy Ek- In 
models with a more general interaction, a k— dependent 
Hartree-Fock shift is also present which, for onsite re- 
pulsion, reduces to a constant Hartree shift mil. The 
selfconsistent fi is the only unknown. 

The high energy expansion about the center of gravity 



G(k,u) = 



M, 



M, 



UJ — U>1 (LO — UJl) 3 (uj — LUl) 4 

Its coefficients 



de ■ (e - wi) A A(k,e); A = 2,3. 



(7) 



(8) 



are called the central moments (Mi = 0, by definition). 
They can be related to correlation functions which occur 
in the short time, or Liouville expansion of the operator 
Cfecr(i) and are evaluated in the limit t = 0. 

It is remarkable that the variance S2 of A(k, e), defined 
by the second central moment 



M, 



s 2 , = m(l — m)U 2 



(9) 



is fc— independent in any dimension d, not only d = oo. 
All the terms in the high energy expansion are sensitive 
to the low energy sector, be it only via the selfconsistent 

We now turn to the CF expansion which is closely re- 
lated to the moment expansion. Formally, it is initiated 
by using u)% and S2 to write the GF as 



G(k,u>) 1 = uj — uji — s\G\(uj). 



(10) 



In this identity, Gi(uj) is again a Herglotz function with 
asymptotics Gi(uj) ~ 1/uj. Iterations, pushing the CF 
further down step by step, require knowledge of the cen- 
ter of gravity OJ21-1 and the variance S21 of Gi-%(cj), to 
write 



Gi^i(uj) 1 =uj - LU21-1 - s 2 21 Gi(uj),1 = 2, 3 . 



(11) 



The two new expansion coefficients depend only on the 
central moments M\ up to the order A = 21 — 1 and 
A = 21 of their respective index. 

By truncating the CF, i.e. by setting Gi(uj) = 0, an 
approximation to the GF is obtained that has / — 1 ze- 
ros and I poles on the real axis. This is defined 2 ^ as 
the PA (I — l\l). It represents the optimal use one can 
make of a set of known spectral moments up to M2i-\- 
The present task, constructing the GF, is rendered essen- 
tially more difficult, because the moments themselves are 
not yet known. A solution based on a PA can be made 
selfconsistent but the moments turn out to be numeri- 
cally quite inexact. This fact is often ignored when it is 
claimed that a certain high energy approximation obeys 
a set of "exact" sumrules. 

We now discuss some well known results concerning 
approximations at the second stage of the CF. As a still 
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exact representation of the GF we have 
G{k,u) = j- 



(12) 



w — lo 3 — s\G 2 (uj) 
The relations between the first few terms are: 
si = M 2 

u 3 = wi + M 3 /M 2 

si = M 4 /M 2 ^M 2 -{M 3 /M 2 f. 



(13) 



Besides the variance, quantities used to further charac- 
terize the internal shape of a spectrum are the skewness 
7 = A/ 3 /m| /2 and the kurtosis k = M 4 /Mf -3. In terms 
of these, we have ui 3 = w\ + js 2 and s| = s 2 (k + 2 — j 2 ). 
From the third moment one finds the coefficient 



lu 3 = (l-m)U + B 3 



H. 



(14) 



This CF coefficient is the first quantity in the expansion 
with a non trivial fc-dependence. The full correlation 
function appearing in the third moment was first derived 
in Ref. 28j and determined selfconsistently for a short lin- 
ear chain in Ref. |29j The shift in the spectral skewness, 
caused by B 3 , regulates the dynamical weight transfer 
between the Hubbard peaks at finite U ^ One can de- 
compose 



B 3 = W + W 3 (k) 



(15) 



in such a way that the term W 3 (k) vanishes in high di- 
mensions. For making contact with the DMFT we will 
presently neglect it and adopt the expression 28 



Bs = 



(2m - 1) 
2m(l — m) 



(16) 



by which it is linked selfconsistently to the expectation 
value of the kinetic energy (T) . 

Concerning the behavior of the fourth moment, not 
even the correlation functions involved in its selfconsis- 
tency loop have as yet been evaluated. Again, the actual 
numerical value of s| is also expected to be sensitive to 
the low energy sector and, in low dimensional systems, 
fc-dependent. 

Given this situation, approximations on the level of 
Eq. (|12f> are at present inevitable. Straightforward trun- 
cation, G 2 (oj) = 0, leads to the PA (1|2). This solu- 
tion with two Dirac peaks goes beyond Hubbard-I, be- 
cause the dynamical weight transfer is taken into ac- 
count. The first example of an approximate resumma- 
tion of the CF is the alloy analogy, developed in the pa- 
per called Hubbard-III. 31 Following Hubbard's notation, 
we approximate s 2 z G 2 (u>) by a /^-independent TF Qh(u>), 
which has to be a Herglotz function. 

The alloy analogy satisfies at least the task (i) of a 
TF, namely to generate finite damping. Far away from 
/i, where the excitations are incoherent, it actually repre- 
sents a physically correct picture. We therefore keep the 



result SIh(uj) — > iD for large u from Hubbard-III. The 
physical reason, why the damping is of the order of the 
bare bandwidth D is that the mean free path is as short 
as one lattice constant. In practice, we incorporate the 
high energy damping in an effective uj 3 



LU 3 — LU 3 + iD 



(17) 



and henceforth deal with a terminator that decays as 
1/uj. This way, we conserve the sumrules, encapsuled in 
the central moments Mo to M 3 . Since Hubbard-III is un- 
realistic at low energies, we do not pursue it any further. 
Nevertheless, it should be noted that Hubbard-III gener- 
ates a branchcut in Q.h{oS), causing the imaginary part to 
drop back to zero and a correlation gap with sharp edges 
to appear, at least in the zero temperature limit. This 
property of Hubbard-III is also not expected to survive in 
improved approximations for the metallic phase. We will 
address the consequences that the absence of a branchcut 
has for the shape of the DOS, both in the CFM and in 
the DMFT. 

To sum up, our approximation to the GF is formally 
similar to Hubbard-III, 



1 



(18) 



u — uj 3 — Q(uj) 



but with a TF, f2#(w) = O(oj) + iD, that retains the 
strong damping of the alloy analogy only at high energy. 
Two successive implementations of the TF with appro- 
priate FL properties at low energy, are the subject of the 
following sections. 



IV. LOW-ENERGY PART 

A FS discontinuity is strictly realized only in the zero 
temperature limit and in a system with no residual dis- 
order. Since T = solutions are hardest to obtain with 
DMFT and, on the contrary, easily implemented with our 
method, we concentrate in the following on this limit. We 
write the standard microscopic definition of a selfenergy 
as a complex correction to the bare excitation rjk'- 



G(k, Cj) 1 = LJ - Tjk - £(fc, U)) 



(19) 



and compare with the inverse of Eq. i|18[l . The high en- 
ergy limit S(fc, oo ) is the difference between two disper- 
sive quantities. In the present case, Eqs. and (JSJ have 
identical dispersion and 



pi = rjk - ui = [i - no 



mil 



(20) 



is, in fact, constant. Within the other approximations, 
discussed in the preceding section, we then obtain the 
fc-independent selfenergy 



E(w) 



-Pi 



(21) 



5 



In this case, as in the DMFT, the FS has the exact shape 
of the uncorrelated system. It is given by all A;— points 
where rjk = in Eq. (JjjJ. The QP peak of weight Z 
at the Fermi level and the step of amplitude Z in the 
momentum distribution are fixed by the conditions 



The solution 



E(0) = 



and 



(0) 



1-1/Z < 0. 



(22) 



(23) 



At finite T or in the presence of a residual diffusive mean 
free path, ImS(O) remains finite. 

Guided by the insight that the strong coupling peak 
is distinct from the Hubbard peaks, we can formulate a 
minimal ansatz for the TF — as 



(24) 



Adding a new stage to the CF is the proper way to " add" 
a pole to the GF. When this TF is inserted in Eq. (fTHft . 
it generates a GF with three complex zeros in the de- 
nominator and two zeros in the nominator, i.e. the same 
structure as the PA (2|3). The connection of the parame- 
ters s\ and <I>5 to central moments M4 and M5 is lost. In 
fact, the very existence of moments beyond M3 has been 
sacrificed by admitting Q3, S4, and uj$ as complex quan- 
tities. They now have to be determined from conditions 
O and |2%|). 

For the Herglotz property one finds 



(Ims4) 2 /Imo;5 < Imcu^ = D 



(25) 



as a necessary and sufficient condition. This causes all 
three poles to lie in the upper half-plane. Further, it 
guarantees a normalized, positive semidefinite A(k, e), 
which also implies quite intricate relations between the 
complex residues. 

Now, the important point is the following: This simple 
ansatz is so heavily constrained by sumrules that it offers 
a selfconsistent solution of the problem, without any free 
parameters. It remains to substantiate this claim and 
then to discuss the quality of the solution. 

After inserting Eq. in Eq. , the conditions lT?2l 
and (|23[) can be brought into a system of two linear equa- 
tions for the unknowns s\ and 0)5. The determinant of 
this system is 



deto 



as\ 



(26) 



and the Herglotz property requires deti > 0. This is a 
constraint on the QP weight Z: In stead of Z < 1 (Pauli 
principle) we have Z < s|/(s2 +P?)- Closer inspection 
reveals that it means the QP cannot take more spectral 
weight than the peak in the PA (1|2) that is nearest to 
p. Since around half filling this weight stays above 1/2, 
it is indeed only a weak constraint. 



4 ~ det 2 



and 



w 5 



P1P2 
det-r 



is expressed in terms of the complex quantity 



Vi 



-(piw 3 



(27) 



(28) 



(29) 



It fulfills the Herglotz condition 1251) with the equality 
sign. This is a consequence of our strong T — con- 
straint E(0) = 0, concerning both the real and imagi- 
nary part. The self energy is now parametrized, up to Z, 
which remains free within a restrained interval and will 
be determined by minimizing the total energy. 

We note, before closing this section, that Ref. Il3l allows 
to define one-pole TF's for the more general case of a 
truncated GF that is expressed as a higher order PA. 
The general algorithm is given, by which the Eqs. 122|) 
and 1|23[) can be fulfilled. 

V. NUMERICAL PROCEDURE 

The uncorrelated chemical potential as function of the 
filling, pa(n), depends only on the kinetic energy part 
and is determined once for all. The DOS per lattice site 
in the U = limit 



p (e) = -ImF (e 
is obtained from the onsite GF 



i0 H 



1 



k 



Mo - Ek 



(30) 



(31) 



A factor two comes from summing over spin directions. 
The DOS of the correlated system 



p{e) = —lmF(e 

IT 



i0 H 



(32) 



is obtained from F(u>) = G(k,w), the on-site GF in real 
space, which is independent of the site index. For a k- 
indepcndent selfenergy such as Eq. 121|) . the on-site GF's 
F(u>) and F (uj) are related to each other by 



F(w) = F (u - E(w)). 



(33) 



The fc-summations can then be carried out by using the 
analytic function that represents the solution for Fq(ui) 
in the limit N — > 00. 

We now turn to the discussion of the selfconsistency 
loops. The condition for p is implemented at T = by 
the integral 



dep(e) 



(34) 
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-0.182 



-0.183 




-0.188 




FIG. 1: Energy minimum E to t(Z) for U = 4, D — 1, and 
p,a = —0.166 corresponding to n — 0.79. The choice of D = 1 
means that all energies are measured in units of the halfwidth 
D of the uncorrelated band. 



According to Eq. i|16|) . the term B$ from the third 
moment requires the selfconsistent determination of the 
kinetic energy, (T) = 2 J_ deA(k,e)E k . One finds 



(T) = — de Im(u)F (uj - Efuj)) - 1) 

* J-oo 

uj = e-iO++Mo-S(e-iO+). 
Finally, the total energy is 



(35) 



Etot = \((T) + 



de e p{e)). 



(36) 



The integrals in (|34|l - (|36(l are carried out numerically. 
For the calculations in this paper we took the on-site GF 



(37) 



In the context of d = oo, it is the GF for a Bethe lattice. 
A halfwidth D = 1 is now used as energy unit. For the 
Herglotz property, it is important to choose the square 
root with a positive real part. The model DOS belonging 
to this GF, 



A) 



(e) = -v/l-(e + /i ) 

7T 



(38) 



is the semi-elliptic function which was also used by Hub- 
bard. 

While searching for the selfconsistent p and (T) at a 
given input n and U , the renormalization Z is still kept 
as a parameter, only limited by the condition deti > 0. 



FIG. 2: Spectral density, comparison of po(e) (dashed line) 
and p(e) (full line) for the same parameters as in Fig. Q 



With these constrained solutions for the GF, we calcu- 
late the total energy E to t- As shown in the example of 
Fig. ^ E to t has a well defined minimum as a function of 
Z. Taking the value which minimizes E to t fixes the last 
parameter Z and defines our solution for the GF. 

The DOS obtained for the same input as in Fig. ^ is 
shown in Fig. [5J together with the U = limit. The 
QP band has the same intensity at the Fermi level as the 
uncorrelated band, p(0) = po(0). This invariance signals 
the unitary limit for the Kondo resonance in the limit 
T = 0. Thus, the reduction of the QP weight does not 
show up in p(0) but in the bandwidth, which is scaled 
down by Z . In a lattice system, this one-to-one relation- 
ship between QP weight and bandwidth only holds when 
the selfenergy is local (/c-independent). 

On the fc-resolved level, near the FS, the QP-pole in the 
complex plane has a parabolic trajectory parametrized by 
Vk, Eq. ©: 



uj*(k) = Zr) k + iT k , 



(39) 



with a scattering rate 



T k = (Z Vk ) 2 /A*. (40) 
The halfwidth for coherent states within the QP band is 



A* 



Zs 2 2 \p 2 \ 2 



D{{l-Z)sl-Zplf 



(41) 



This formula is well behaved also in the weak coupling 
limit, in fact for all possible metallic, unpolarized regimes 
of the Hubbard model. In the strongly correlated regime, 
the energy scale A* is smaller than ZD, so that the 
excitations in the wings of the QP band cease to be 
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coherent. Since we have modeled the ballistic limit 
(residual diffusive scattering rate = 0) the QP res- 
onance in A(k, e) is a Dirac peak for k = kp and, for 
k kF , it has the so called Breit-Wigner lineshape (see 
Hcdin and Lundquist- 2 for a generic plot). This shape is 
due to an interference between the QP residue and the 
other residues. The lineshape becomes approximately 
Lorentzian whenever a T ( i > Tk is present. 

Returning to the DOS, we note that the global shape of 
the valence spectrum for a hole doped Mott insulator, i.e. 
QP band and lower Hubbard band, is well rendered by 
our present approximation. The sumrules up to M 3 are 
exactly satisfied and their interplay regulates the overall 
skewness and the relative weight of all three features. 

The one-pole TF has the drawback of being unable to 
reproduce a sharp gap formation. The high level of in- 
tensity between the QP band and the upper Hubbard 
band shows that the dynamical spectral transfer— is not 
realized completely, at least for U/D = 4. The intensity 
at the minimum decays like (U/D)~ 2 , so that this spu- 
rious effect disappears for larger U. We shall discuss the 
presence of residual intensity in the gap region in more 
detail when we compare with our second ansatz and with 
the DMFT. 

Two remarks to conclude this section: (i) The limit 
U = oo describes spin- and charge-excitations in the sub- 
space of singly occupied sites. It is equivalent to the 
t — J- model with J = 0. The one-pole TF thus allows 
to project out a quantitatively valid GF for the valence 
sector near this limit, up to terms of order (U/D)~ 2 . 
(ii) Ratios U/D rs 4 are relevant for the doping controlled 
Mott transition in real materials close to criticality. Our 
main motivation to pursue the CFM was to investigate 
whether by simply adding a second complex stage to the 
TF we could handle this regime in a semiquantitative 
way. The derivation of the two-pole TF and its applica- 
tion to U/D = 4 are presented in the next section. 



VI. IMPROVING THE DYNAMICAL WEIGHT 
TRANSFER 

To generalize our ansatz, we introduce algebraic ex- 
pressions for f2(w), such that Eq. (|18|l can be cast into 
the form of a truncated CF with complex coefficients. 
This defines the CFM, provided the Herglotz condition 
is satisfied. The fc-resolved GF has then the structure 
of a generalized higher order PA. By terminating the PA 
(1|2), we still retain the important sumrules that govern 
the dynamical weight transfer. Spectral moments beyond 
M3 cannot be recovered, but this may not be a great 
sacrifice, given the difficulties known from the projection 
method to obtain correct values for higher moments. 

What can be gained by using complex coefficients is the 
possibility to model constructive and destructive interfer- 
ence phenomena in the GF at intermediate energies. A 
single feature in the spectral function can be built up by 
the contributions of several poles, resulting in uncommon 



lineshapes. An ansatz frequently employed in the phe- 
nomenological interpretation of spectra is the superpo- 
sition of complex poles with real residues (superposition 
of Lorentzians in the spectrum). Although this allows 
several peaks to coalesce, it still eliminates interference. 
One striking example of a Fano like interference within 
the coherence range of halfwidth A* is the Breit-Wigner 
lineshape of the QP— As discussed in the preceding sec- 
tion, the complex residues of the two valence poles in the 
GF (arising from the one-pole TF in the large U limit) 
are enough to obtain this lineshape. 

Likewise, the dynamical weight transfer and the for- 
mation of the correlation gap can be interpreted as a 
destructive interference in the intermediate energy range 
between the Hubbard bands. When the interference is 
complete the function G2(oj) in Eq. l|12f> should acquire 
a branchcut and a gap interval with zero DOS and sharp 
edges should result. This may be possible only on the in- 
sulating side of the Mott transition and strictly at T = 0. 
When the system is metallic and the chemical potential 
falls in a region of high DOS, it is satisfactory to model 
the correlation gap by a deep minimum. We demonstrate 
here that this situation is captured by a two-pole TF of 
the form 

fl(w) = 5— . (42) 

U) — LO5 

UJ — Ul'J 

The new degrees of freedom are given by s\ and 0)5. 
These will be found due to some qualitative arguments, 
restricting the ansatz from the start. Then, s| and w-j 
can again be eliminated by the FL conditions of Eqs. 
(|22|l and (12;it , using the next iteration of the algorithm 
in Ref. El 

The GF now has four poles and the Herglotz condition 
becomes a crucially important issue. To formulate it, for 
arbitrary complex values of <I>3 to ujy , seems at first sight 
rather difficult. The GF on the FS (Eq. (JTHJ) with rjk = 0) 
has additive coherent and incoherent contributions, 

1 =Z +Gb{uj) , (43 ) 

u) — n(u>) ui 

The decomposition is possible, because one pole lies ex- 
actly on the real axis. This will enable us to manage the 
Herglotz condition for £(o>) more easily: from Eq. 142|) 
we obtain a background function Gj,(lo) with three poles 
that can be written 

G 6 M = J—^-g . (44) 



The new coefficients are designated by capital Greek let- 
ters. Systematically, they depend on Z and, at order A, 
on all coefficients in Eqs. (|18|l and (|42|l with index A < A. 
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Explicitly, the first three are 
Pi 



fll = - 



1 - Z 
Z det 2 

1 \ a p 2 s\ 
Pi 



(45) 



dety 



deto 



in terms of the previously defined quantities, Eqs. 10, 
(EU, (PI), and @. 

The high energy damping in the background function 

is 



(1 



det-. 



)Imw3 > Imw3 = D. 



(46) 



Between the coherence energy A* in Eq. (|4U|) and the 
background function at the Fermi edge there is the rela- 
tion 



A*ImG 6 (0) = Z. 



(47) 



For £4 = £^5 = 0, we recover the one-pole TF and 
Eq. 141|l for A*. Since f2i and S| are reai J there are now 
only three complex quantities and the Herglotz condition 
can be specified exhaustively, analogous to Eq. (|25|l : 



(ImE^/Im^ < Imft 3 . 



(48) 



The foregoing analysis suggests that £4 and are more 
useful than s 4 and 0)5 as control parameters. To obtain 
the selfenergy, one can then use Eqs. i|43[) and (|44|) . For 
further discussion we parametrize 



ft 3 = X 3 + iY 3 
£4 = X 4 + iY 4 
^5 = X 5 +iY 5 - 



(49) 



A minimal requirement for properly defining the dynam- 
ical weight transfer is vanishing Im£(e — i0 + ) in one other 
point e = Xq on the real axis, apart from e = 0. It hap- 
pens if (and only if) the equality sign applies in iffify . 

Y%Y 5 and to 



This leads to the condition Y 4 



x = X 5 



Y5 
Ya 



-Xa 



(50) 



for the position. The Herglotz property guarantees that 
it is in fact a minimum. The influence of this interference 
on the shape of the valence spectrum is weakest for X$ — 
0. For simplicity, we also need to set Y5 = Y4 = Y 3 , where 
Y3 is already defined in Eq. (|46|l . The last parameter 
X 4 = — xo is then fixed by the point with lowest intensity 
inside the correlation gap. 

Before continuing with this ansatz, it is important to 
realize that it cannot apply exactly at half filling. There, 
the metallic phase is obtained by driving U/D below the 
critical ratio (so called bandwidth controlled transition) X 
The particle-hole symmetric DOS has a quite different 




FIG. 3: Spectral weight Z of QP pole as function of electron 
density n for U = 4. Comparison of the Gutzwiller approx- 
imation (GA) to our result with the one- and two-pole TF's 
(CFM1 and CFM2). 



morphology than what is shown in Fig. [3 the QP's are 
in the center and the correlation gap is split in two sym- 
metric gaps of order U/2£ In our approach, it can be 
envisaged to use one additional complex stage to model 
two symmetric destructive interferences. 

In the doping controlled regime, there is only one large 
correlation gap and we have a good qualitative argument 
for xq: the strong skewness (large |pi|) causes the QP 
band and the minimum position to always be on op- 
posite sides of the center of gravity. This is well satisfied 
by setting 



.To = Recc>3. 



(51) 



The remaining free parameter Z is determined again by 
minimizing the total energy. The numerical procedure is 
as described before. In Fig. |3 results with the one- and 
two-pole TF's (CFM1 and CFM2) are compared to the 
Gutzwiller approximation (GA) at constant U , as func- 
tion of the filling. The upper curve is the well known 
lower bound for the GA, Z = (1 — n)/(l — to), obtained 
by excluding double occupancy. By projecting out the 
background, the GA is known to systematically overesti- 
mate the coherent weight. The behavior that results from 
the CFM, i.e. : lowering of Z and upward curvature at 
the approach of zero doping (1 — n — > 0), is close to that 
of the exact Kondo scale in the Bethe ansatz solution for 
the Anderson impurity. 19 



VII. COMPARISON WITH OTHER METHODS 

With the two-pole TF, realistic results for the DOS 
in the doping controlled regime can be obtained, even 
close to the critical U. To illustrate this, we compare our 
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CFM with the DMFT for two different impurity solvers. 
The impurity solvers perform the crucial step in mapping 
the Hubbard lattice model onto an Anderson impurity 
model. The effective medium surrounding a given site 
is determined self-consistently, still a formidable many- 
body problem. The NRGj^S used to solve it at the lowest 
temperatures and energies, requires a heavy amount of 
computer time. The NC A 21 i 22 is an alternative, more 
analytic method, less reliable for |w| <C A*, but obeying 
high energy sumrules well. Therefore, NRG and NCA 
are expected to be complementary. 

A comparison for the same parameters as before, i.e. 
U = 4 and n = 0.79JS shown in Fig. H The NRG data 
are taken from Ref . fla. NCA is our own unpublished cal- 
culation, CFM1 is again the DOS from Fig.^land CFM2 
the result with Eq. I|42[l . All four solutions obey the 
p(0) = po(0) condition. This confirms that temperatures 
in the DMFT solutions are sufficiently low to warrant 
a comparison with our T = results. As manifest in 
the width of the QP band, the selfconsistent Z obtained 
for CFM2 coincides with both versions of DMFT. Since 
NRG is expected to determine essentially the exact low 
energy scale, this is a good point for both the NCA and 
the CFM2 results. 



/i and i?3. The quantitative agreement with the DMFT 
in the QP band and good overall agreement in the entire 
valence sector is due to this built in interference. In com- 
paring CFM1 and CFM2, one notices a feedback of the 
improved gap region on the QP band: The sumrules up 
to M3 are satisfied for both approximations, but the dy- 
namical weight transfer is more complete within CFM2. 
Removing the spurious intensity inside the gap slightly 
raises the QP weight (Compare Fig. I2J), bringing it in 
agreement with the NRG. 

The rather large variation among the different solu- 
tions in the region of the upper Hubbard band is re- 
markable and still deserves more detailed investigations. 
At higher temperatures, T > A*, where Quantum Monte 
Carlo is available as benchmark, the NCA was found to 
be satisfactory^ In the present comparison, the NCA 
comes closer to obeying the sumrules than the NRG. As 
far as numerical effort is concerned, the NRG is the most 
demanding, followed by the NCA. The CFM2 stands up 
quite honorably in this comparison, especially when con- 
sidering that the sumrules are rigorously incorporated, 
no "technical" broadening needs to be introduced and 
the required computer time to achieve selfconsistency is 
in fact negligible. 




Energy 

FIG. 4: Spectral density. Comparison of the continued frac- 
tion method using the one- and two pole TF's (CFM1 and 
CFM2) with NCA and NRG (data taken from Ref. [3). 

The solutions start to differ somewhat in the gap re- 
gion. Neither DMFT version shows a gap with sharp 
edges that would correspond to a branchcut in the self- 
energy. A real benchmark for low T impurity solvers in 
the doping controlled regime does not yet exist. From 
the NCA, we can confirm that some very low residual 
density inside the gap seems to be the generic situation. 

In the ansatz for CFM2, the existence of a point with 
zero DOS is postulated. Determining its position accord- 
ing to Eq. H51f) involves the selfconsistency conditions for 



VIII. DISCUSSION AND OUTLOOK 

We present the CFM as a new method to calculate the 
selfenergy, as well as various k— resolved and (partially 
or fully) k— integrated spectral functions of the Hubbard 
model in the correlated metal phase. We expand the sin- 
gle particle Green function as a continued fraction, as far 
as moment sumrules are exploitable, and then use a prop- 
erly chosen terminator function. In this paper, moment 
sumrules up to M3 are implemented and the "termina- 
tor" is a fc— independent complex function with one or 
two poles that obeys the correct Fermi liquid properties 
at low energies. In this local approximation to the selfen- 
ergy, we compare our results for the density of states with 
the DMFT. Our method has a precision comparable to 
state-of-the-art impurity solvers NRG and NCA. It cov- 
ers all energy scales reliably, whereas the low T impurity 
solvers each have their strengths and weaknesses. 

With the second stage in the terminating function we 
are able to improve the dynamical weight transfer be- 
tween the upper and lower Hubbard peak and thereby ob- 
tain very good agreement with DMFT for the QP weight 
Z or low energy scale. This is significative, because NRG- 
DMFT yields the exact result for this quantity. Unlike 
the time consuming DMFT calculations, the CFM uses 
simple, algebraic functions, for which selfconsistency con- 
ditions are rapidly found. 

The CFM is generalizable in many directions. How- 
ever, the possibility to circumvent heavy manybody cal- 
culations by such a simplified ansatz seems too attractive 
to be true. Thus, before advocating possible extensions, 
we need to analyze the reason for the quantitative success 
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of the CFM in the strong coupling limit. The Hubbard 
model with a local selfenergy is, admittedly, only a toy 
model but nevertheless an obligatory testing ground for 
this important issue. 

The algebraic terminator functions were already in- 
troduced earlier. Their Fermi liquid properties, essential 
for circumventing the explicit manybody calculations, are 
determined by using the Luttinger sumrule as an input. 
Their phenomenological possibilities could be demon- 
strated by leaving Z as a free parameter li&iLiSi The cor- 
nerstone of the CFM as a microscopic method is now the 
variation of the total energy to obtain Z . Given G(k, ui), 
we calculate the total energy from the exact manybody 
expression, actually another sumrule first found by Gal- 
itski. However, without an explicit wavefunction, we 
have no rigorous variational principle. In making the 
Gutzwiller approximation, beyond the Gutzwiller ansatz 
for the wavefunction, one is also abandoning the rigor- 
ous variational principle but one keeps Z as variational 
parameter. 

The answer to the question, why our method is vari- 
ational, is probably that we are using a GF, fully con- 
strained by sumrules, that leaves no other free parameter 
but Z. To obtain the Kondo effect, we need degeneracy. 
Our model has spin degeneracy, Nf = 2, for its flavors. A 
clue, why including the incoherent background spectrum, 
instead of projecting it out, improves the outcome for Z 
comes from the limit of large Nf&m The low energy scale 
(Kondo temperature) in the Anderson impurity model 
can be obtained exactly, as function of n and m = n/Nf. 
Also, coherent spectral weight is of order zero in 1/Nf, 
the leading background contribution starts at first order. 
Neglecting background, as for instance in the slave boson 
method at mean-field level, yields Z = (1 — n)/(l — to), 
as plotted for Nf = 2 in Fig. [3] Compared with the 
Bethe ansatz, this renormalization is not enough. Now, 
the influence of the background is strikingly illustrated 
by solving for the Kondo temperature only to the first or- 
der in This causes indeed a substantial decrease, 
bringing the result close to the exact value. The cor- 
rect doping dependence displays the upward curvature, 
as also seen in our approximations CFM1 and CFM2. 
Finally, the improvement from CFM1 to CFM2 shows a 
delicate interplay between the dynamical weight trans- 
fer, related to the double occupancy, and the low energy 
scale. 

If determining Z by varying the total energy is indeed 
a valid variational principle, it makes the CFM indepen- 
dent of the limit d — oo, thus giving it high flexibility and 



a large field of applications. It is straightforward and, for 
low dimensional systems, potentially very important to 
incorporate the k— dependence in the moment M3. The 
term W^,(k) in Eq. (|15f) was already identified in the ex- 
act diagonalization of a short linear chain^S as causing 
a coupling of the QP to antiferromagnetic fluctuations. 
This can be generalized to fluctuations above other pos- 
sible groundstates and the selfconsistent determination 
of Wa(k) thus offers a path to describing the feedback of 
bosonic fluctuations on the low energy sector. Up to now, 
the treatment of low energy effects within the projection 
method was based more on physical intuition, or guess- 
work for the more critical observer, than on an objective 
procedure. 

The extension of the CFM to higher moments becomes 
possible due to its close relationship with the numerical 
Lanczos procedure for finite lattices. All what is missing 
is a proper termination of the continued fraction with a 
TF representing the low energy sector and the dissipa- 
tion. The general algorithm for calculating the coeffi- 
cients in the TF is given in Ref. PHI . 

As an outlook, we enumerate other possibilities that 
are inherent in the CFM, beyond the results of this paper. 
They are listed roughly according to increasing effort that 
will be required to implement them, (i) A more detailed 
exploitation of spectral functions on the /c-resolved level: 
e.g. the interpretation of Raman, ARPES, or tunnel- 
ing data requires the partial summation of A(k, e) over 
selected spots in the Brillouin zone, weighted by ma- 
trix elements, (ii) Hubbard lattice models with a more 
realistic kinetic energy part, including van Hove singu- 
larities, (iii) The generalized periodic Anderson model 
(PAM): lattice models with Hubbard repulsion among 
transition orbitals but, in addition, hybridization with 
ligand orbitals. (iv) Implementation of LDA+CFM. The 
algebraic simplicity of the CFM allows to calculate the 
charge transfer effects, present in model Hamiltonians of 
the PAM type, on an "ab-initio" level. These effects, im- 
portant for many real materials, could not yet be handled 
successfully by LDA+DMFT. (v) Not difficult to imple- 
ment, but leaving the strict framework of the CFM as 
an algebraic method, is the inclusion of non Fermi liquid 
effects on a phenomenological levels 

In conclusion, we have attempted to demonstrate by 
means of the Hubbard model that the CFM is a powerful 
method. Numerically simple, due to its algebraic struc- 
ture, it is still sufficiently rigorous to deal with strongly 
correlated electrons in mesoscopic and macroscopic sam- 
ples of condensed matter. 
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